library(sdmTMB)library(tidyverse)library(tidyverse)library(viridis)library(devtools)library(terra)library(tidylog)library(RColorBrewer)library(patchwork)library(egg)library(ggpmisc)library(ggh4x)library(ggtext)# Set pathhome <- here::here()# Source functions for overlap, predation and plottingfor(fun inlist.files(paste0(home, "/R/functions"))){source(paste(home, "R/functions", fun, sep ="/"))}# Palette for plottingpal <-brewer.pal(name ="Paired", n =8)
Load prediction grid
# scale environmental variables with respect to catch data, and diet covariates with respect to diet datacatch <-read_csv(paste0(home, "/data/clean/catch_clean.csv")) |>drop_na(depth, oxy, salinity, temp)diet <-read_csv(paste0(home, "/data/clean/stomachs.csv"))pred_grid <-bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv"))) |>filter(quarter ==4) |># Not needed in theory for saduria...drop_na(saduria) |>mutate(year_f =as.factor(year),quarter_f =as.factor(quarter),month_f =as.factor(month),oxy_sc = (oxy -mean(catch$oxy)) /sd(catch$oxy),temp_sc = (temp -mean(catch$temp)) /sd(catch$temp),temp_sq = temp_sc^2,depth_sc = (depth -mean(catch$depth)) /sd(catch$depth),depth_sq = depth_sc^2,salinity_sc = (salinity -mean(catch$salinity)) /sd(catch$salinity),pred_length_sc =0,saduria_sc = (saduria -mean(diet$saduria))/sd(diet$saduria))
Predict with sims from density and diet models to propagate uncertainty
# predict and then for loop to summarise and avoid having too large pred grids in the memorynsim <-500sim_dens <-predict(mcod, newdata = pred_grid, nsim = nsim) |>as.data.frame()sim_her <-predict(mher, newdata = pred_grid, nsim = nsim) |>as.data.frame()sim_spr <-predict(mspr, newdata = pred_grid, nsim = nsim) |>as.data.frame()sim_sad <-predict(msad, newdata = pred_grid, nsim = nsim) |>as.data.frame()
The density prediction is used later for spatial overlap + it gets left_joined to the other sims, so it gets a special treatment
# A prediction grid cell is 3*3 kmarea <-9# Summarise predator biomass by year and sim and join that to calculate per capita predationpred_dens <- pred_grid_density |>summarise(cod_biomass =sum(sim_density*area),.by =c(year, sim)) # Apply functions and combinetot_pred <-bind_rows(get_annual_predation(sim_her, threshold =1, power =2, prey_name ="Herring"),get_annual_predation(sim_sad, threshold =1, power =2, prey_name ="Saduria"),get_annual_predation(sim_spr, threshold =1, power =3, prey_name ="Sprat"))# Prepare data for fitting gams and for plottingclean_pred <- tot_pred |>pivot_longer(c(-year, -species)) |>mutate(var =ifelse(grepl("_sd", name), "sd", "mean"),metric =ifelse(grepl("cap", name), "Per capita (kg/kg)", "Total (tonnes)")) |> dplyr::select(-name) |>pivot_wider(values_from = value, names_from = var) |>mutate(species =as.factor(species))# Fit gamma-GAMs to annual predation metricsgamma_gam_tot <-sdmTMB(mean ~ species +s(year, by = species),spatial ="off", family =Gamma(link ="log"),data = clean_pred |>filter(metric =="Total (tonnes)") |>mutate(mean = mean/1000))gamma_gam_cap <-sdmTMB(mean ~ species +s(year, by = species),spatial ="off", family =Gamma(link ="log"),data = clean_pred |>filter(metric =="Per capita (kg/kg)"))nd <-expand_grid(species =as.factor(c("Sprat", "Herring", "Saduria")),year =unique(clean_pred$year))p_tot <-predict(gamma_gam_tot, newdata = nd, se_fit =TRUE)p_cap <-predict(gamma_gam_cap, newdata = nd, se_fit =TRUE)p <-bind_rows(p_tot |>mutate(metric ="Total (tonnes)"), p_cap |>mutate(metric ="Per capita (kg/kg)"))clean_pred <- clean_pred |>mutate(lwr = (mean - sd),upr = (mean + sd),lwr =ifelse(lwr <0, 0, lwr),mean =ifelse(metric =="Total (tonnes)", mean /1000, mean),lwr =ifelse(metric =="Total (tonnes)", lwr /1000, lwr),upr =ifelse(metric =="Total (tonnes)", upr /1000, upr))p0 <-ggplot(clean_pred, aes(year, mean)) +geom_errorbar(aes(ymin = lwr, ymax = upr), width =0, color ="grey40", alpha =0.7,linewidth =0.4) +geom_point(color ="grey30", alpha =0.8) +geom_line(data = p, aes(year, exp(est)), color = pal[2], linewidth =0.8) +geom_ribbon(data = p, aes(x = year,ymin =exp(est -1.96*est_se),ymax =exp(est +1.96*est_se)),fill = pal[1], alpha =0.3, inherit.aes =FALSE) +facet_grid2(metric~species,scales ="free_y", independent ="y") +theme(aspect.ratio =6/7,strip.text.x.top =element_text(size =10, margin =unit(rep(0.06, 4), "cm")),strip.text.y.right =element_text(size =10, margin =unit(rep(0.06, 4), "cm"))) +labs(x ="Year", y ="Predation")tag_facet(p0, fontface =1, col ="gray30", size =3)
ggsave(paste0(home, "/figures/pred_yearly.pdf"), width =17, height =9, units ="cm")
Plot weighted predation in space
spatial_pred <-bind_rows(get_spatial_predation(sim_her, years =c(1994, 2022), threshold =1, power =2, prey_name ="Herring"),get_spatial_predation(sim_sad, years =c(1994, 2022), threshold =1, power =2, prey_name ="Saduria"),get_spatial_predation(sim_spr, years =c(1994, 2022), threshold =1, power =3, prey_name ="Sprat"))
pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]
mutate: new variable 'annotateText' (character) with 2 unique values and 0% NA
p1 + p2 + p3 +plot_layout(axes ="collect")
ggsave(paste0(home, "/figures/spatial_predation.pdf"), width =19, height =13, units ="cm")
Plot CV in space for cod density and predation
# Apply functions and combinecv_pred <-bind_rows(get_cv_predation(sim_her, years =c(1994, 2022), threshold =1, power =2, prey_name ="Herring"),get_cv_predation(sim_sad, years =c(1994, 2022), threshold =1, power =2, prey_name ="Saduria"),get_cv_predation(sim_spr, years =c(1994, 2022), threshold =1, power =3, prey_name ="Sprat"))
pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]